clc; clear;

% 参数范围和采样数
x110_list = linspace(-0.7, 0.2, 150);  % x110参数扫描
q110_0 = 1;      % 固定初始条件 q110
q210_0 = 0;      % 固定初始条件 q210
k = 0.2;         % 固定耦合强度
phi0 = 0;        % 固定初始phi

% 其他固定初始条件（x120,x130,x210,x220,x230）
x_init_base = [0; 0; 0; 0; 0];

dt = 0.01;
T_total = 600;
T_trans = 550;
steps = round(T_total/dt);
trans_steps = round(T_trans/dt);

% 参数矩阵定义（参考你代码）
alpha = 2.4;
A = [-alpha, 0, 0;
      0, 0, 0;
      0, 0, 0];
B = [1, -4, -3.5;
     0, 1, 2;
    -1, -4, 1.5];

% 结果存储
q11_peaks_all = cell(length(x110_list),1);
phi_peaks_all = cell(length(x110_list),1);

parfor i = 1:length(x110_list)
    x110 = x110_list(i);
    % 构造完整初始状态向量 [x110,x120,x130,x210,x220,x230,q110,q210,phi]
    x0 = [x110; x_init_base; q110_0; q210_0; phi0];
    
    % 定义ODE函数句柄，调用你的mCNN函数格式
    f = @(t,x) mCNN(t, x, A, B, k);
    
    % RK4积分
    X = zeros(length(x0), steps);
    X(:,1) = x0;
    for n = 1:steps-1
        k1 = f(0, X(:,n));
        k2 = f(0, X(:,n) + dt/2*k1);
        k3 = f(0, X(:,n) + dt/2*k2);
        k4 = f(0, X(:,n) + dt*k3);
        X(:,n+1) = X(:,n) + dt/6*(k1 + 2*k2 + 2*k3 + k4);
    end
    
    % 丢弃暂态，截取稳态数据
    X_steady = X(:, trans_steps:end);
    
    % 提取q11和phi的峰值
    q11_signal = X_steady(7, :);  % q11为第7维
    phi_signal = X_steady(9, :);  % phi为第9维
    
    q11_peaks_all{i} = peakdetect(q11_signal);
    phi_peaks_all{i} = peakdetect(phi_signal);
end

% 整理绘图数据
xdata_q11 = [];
ydata_q11 = [];
for i = 1:length(x110_list)
    xdata_q11 = [xdata_q11, repmat(x110_list(i), 1, length(q11_peaks_all{i}))];
    ydata_q11 = [ydata_q11, q11_peaks_all{i}];
end

xdata_phi = [];
ydata_phi = [];
for i = 1:length(x110_list)
    xdata_phi = [xdata_phi, repmat(x110_list(i), 1, length(phi_peaks_all{i}))];
    ydata_phi = [ydata_phi, phi_peaks_all{i}];
end

% 绘图
figure;
plot(xdata_q11, ydata_q11, '.r');
hold on;
plot(xdata_phi, ydata_phi, '.b');
xlabel('x_{110}');
ylabel('q_{11},\phi');
grid on;


% 你的mCNN函数，保持不变
function dxdt = mCNN(~, x, A, B, k)
    X1 = x(1:3);
    X2 = x(4:6);
    q11 = x(7);
    q21 = x(8);
    phi = x(9);

    M1 = exp(-q11);
    M2 = exp(-q21);
    Y1 = M1 .* X1;
    Y2 = M2 .* X2;

    W_phi = 1 - 2*exp(phi);
    K = diag([0, 0, k]);

    dX1 = -X1 + A*Y1 + B*X1 + K*W_phi*(X1 - X2);
    dX2 = -X2 + A*Y2 + B*X2 - K*W_phi*(X1 - X2);
    dq11 = X1(1);
    dq21 = X2(1);
    dphi = X1(3) - X2(3);

    dxdt = [dX1; dX2; dq11; dq21; dphi];
end

% 简单峰值检测函数
function peaks = peakdetect(signal)
    [pks, ~] = findpeaks(signal);
    peaks = unique(pks);
    if isempty(peaks)
        peaks = max(signal);
    end
end
